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1 Introduction 

With  the  capacity  of  today’s  computers,  one  can 
envisage  the  resolution  of  shape  optimization  prob- 
lems in  aerodynamics.  Nevertheless,  optimization 
methods  require  many  evaluations  of  different  aero- 
dynamic configurations,  and  so  are  much  more  ex- 
pensive than  a single  analysis.  It  is  therefore  manda- 
tory to  find  methods  that  evaluate  aerodynamic 
functions  and  their  gradient  at  the  lowest  possible 
computational  cost,  as  well  as  fast  and  robust  op- 
timization methods. 

Classical  optimization  techniques  (descent  methods) 
not  only  require  the  value  of  the  function  to  opti- 
mize, but  also  of  its  gradient.  The  classical  way 
to  compute  the  gradient  is  to  use  a finite-difference 
formula;  the  main  drawback  of  this  method  is  due 
to  the  fact  that  n + 1 evaluations  of  aerodynamic 
functions  are  necessary  at  each  iteration,  n being 
the  number  of  parameters  defining  the  geometry 
to  optimize.  So,  such  methods  are  completely  un- 
suited to  aerodynamic  shape  optimization,  because 
of  the  high  computational  cost  of  the  single  analy- 
sis. Alternative  methods  (stochastic  optimization, 
genetic  algorithms)  that  don’t  require  gradient  in- 
formation are  also  highly  costly  in  term  of  CPU 
time. 

For  a few  years,  techniques  for  sensitivity  analysis 
based  on  the  optimal  control  theory  have  been  de- 
veloped ([l]j[13]).  These  techniques  derive  from  the 
state  equations  another  set  of  equations  called  ’’ad- 
joint” or  ’’costate”  equations.  The  solution  of  these 
adjoint  equations  is  used  to  compute  the  gradient 
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at  very  low  cost;  since  solving  the  adjoint  equa- 
tions is  equivalent  to  solve  the  state  equations,  the 
cost  of  sensitivity  analysis  is  greatly  reduced.  Some 
authors  use  adjoint  equations  derived  from  the  dis- 
cretized Euler  equations  ([5], [7]).  In  this  paper,  we 
focus  on  adjoint  equations  derived  from  analytical 
state  equations. 


2 Optimal  control 

A typical  shape  optimization  problem  can  be  stated 
as  follows: 

Let  Q,  be  a subspace  of  R2.  Find  the  shape  of  T o, 
a boundary  of  Q controlled  by  design  variables  u, 
such  that  the  functional : 

J(  u)=  / ^{s,qx,qy,u)dTJ  (1) 

Jr, 

computed  on  a boundary  Tj  (where  T j and  To  can 
be  different)  is  minimized 

Over  the  domain  Q,  the  state  variables  s are  gov- 
erned by  the  advective-diffusive  state  equations: 

dx  [f (s)  - f„(q,  q®,  qy)]  + dy  [g(s)  - g„(q,  q®,  qy)]  = 0 

G,-(q,  q®,  qj/,  u)  = 0 

where  q = q(s)  are  the  primitive  variables  and 

<lx  = dxq=Mdxs  (4) 

q y = dyq=MdyS  (5) 

with  M = ^ the  jacobian  of  the  transformation 

between  primitive  and  state  variables. 


on  Q (2) 
on  Ti  (3) 
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In  order  to  find  an  optimality  condition,  one  have 
to  form  the  Lagrangian  of  the  problem: 


C = J 

+ [ 4>T[dx(f-fv)  + dy(g-g „)]<K1 


+ / Aj  (Qr  - Mdxs)  + Ay  (qr  - M^s)  dQ 

Ja  (6) 


The  local  variation  due  to  the  design  variables  is: 
= (8) 


where  the  Hamiltonians  % are  defined  as: 


• l~tf  — $ + (j/j  Gj  on  r j?  ; 

• Hi  = fij Gt-  on  the  other  boundaries. 


where  ip,  Hi,tpx,ipy  are  the  costate  or  adjoint  vari- 
ables  associated  with  the  state  equations,  state  bound- 
ary conditions  and  primitive-conservative  realtions 
respectively. 

The  optimal  control  theory  states  that  the  Lagran- 
gian  is  stationary  at  the  optimum  solution,  i.e.  at 
the  solution,  the  (independant)  variations  of  the  La- 
grangian with  respect  to  the  state  (s),  costate  [ip,  ji) 
and  design  (u)  variables  is  zero. 

The  variation  of  the  Lagrangian  due  to  the  design 
variables  is  the  sum  of  two  contributions: 

• a local  variation  due  to  the  variation  of  the 
design  variables; 

♦ a convective  variation  due  to  the  displacement 
of  the  domain  and  boundaries. 

The  variations  of  the  Lagrangian  due  to  the  state, 
and  costate  variables  are  local  variations. 


The  local  variation  due  to  the  state  variable  is, using 
Gauss’  theorem: 


&C*  = 


f {ipT  [(A  - A v)nx  + (B  - B„)  ny]}  SsdTi 

Jl ir 

- [ { [\*nx  + A Tyny]  M + pf  dsG{  + ds$]  6s  dTi 

J ur 

+ J [{A-Av)Tdxi>  + {B-Bv)Tdyip}Ssdn 

~i{ 


dx  {MT Xx)  - dy  (MTAy)  + A l —4s 


+A  y^dyslssdQ 


(9) 


SC  ? = 


[ {'ipT  [(A  - Av)  nx  + (B  - Bv)  ny]}  SsdTi 
J ur 

- [ { [A*  n*  + A Tyny\  M + $ dsG{  + <9S$}  6s  dTi 
J ur 

+ J {( A-Av)Tdxip  + {B-Bvfdyil>}SsdQ. 

- f {MTdxXx-MTdy\y}8sdQ  (10) 

Jn 


2.1  Local  variations 


The  local  variation  of  the  lagrangian  can  be  split  in 
eight  parts: 

8C^  = 8C^  + 5C^+5C^+8Cl<  + 8Ct 
+ E^‘+<^  + ^ (7) 

i 

As  the  variations  of  the  state,  costate  and  design 
variables  are  independent,  the  variation  of  their  con- 
tribution to  the  variation  of  the  Lagrangian  must 
be  zeroed  independently. 

It  is  easy  to  demonstrate  that  setting  SCfoc  and  SC £*c 
to  zero  leads  to  the  state  equations  and  bound- 
ary conditions,  while  setting  SC **  and  SCt> yc  to  zero 
leads  to  the  primitive-conservative  relations. 


where 


A 


A y 


df_ 

ds 

dU 


B 


% 

3s 


- 3q 


(ii) 


The  local  variation  due  to  and  qy  are,  using 
Gauss’  theorem: 


+ E $ Gi  “ ( Cn*  + Eny) 

i 


+ 1 {\x+(CTdx^  + ETdyi>)}5qxdQ 

Jn 


SqxdT 

(12) 


6C?Z  = 


where 


-ipT  (Dnx  + Fny)}  dqy  dT 


C = 


of  a functional  of  the  form: 

J = JrfdT 


* + (DT' 

dxrp  + F Tdyip) } £qy  dO 

is: 

with: 

L E 

- 

dfv_  F 

= 

where: 

Sqj, 

d<ly 

(13]jConv=  j ^V/-W  + /c/)  dr 


daCo  * t 

W~ 


Equaling  the  integrals  over  Q to  zero  in  (12)  and 
(13)  gives  the  relations  between  rp)\x  and  Ay: 

A*  = - (CTdxip  + ETdy)  (14) 

Aj,  = ~{DTd^  + FTdy)  (15) 


• the  boundary  T is  described  by  the  parameter 
a:  T : (ar(cr),  J/(cr»; 

• t is  the  tangent  vector:  t = ( dax  day ); 

• Co  is  the  boundary  displacement  due  to  the 
change  of  the  design  variables. 


Equaling  the  integral  over  Q to  zero  in  (10)  and  us- 
ing the  two  above  relations  gives  the  adjoint  equa- 
tions: 

(A  - A vf  dxrj)  + (B  - B „f  dy4> 

+Mt  (dxfv  + dygSj  = 0 (16) 

with: 

f„  = (CTdxip  + ET  dy  ip) 
g„  = (Dt5iV-  + FTdyi>) 
the  ” adjoint  viscous  fluxes”. 

The  adjoint  boundary  conditions  are  found  by  equal- 
ing integrals  over  I\*  to  zero  in  (10), (12)  and  (13): 

(An  - A„„)T  + CTndx^  + E Tndyi> 

+nJdsGi+ds$  = Q (17) 

d^  + tiJd^Gi  = CU  (18) 

dqy$  + dqyGi  = (19) 

with: 

A xi  — A fix  Bny  A un  ~ A.yTix  -h  B^ny 
Cn  = C nx  T Dny  En  = E nx  + F ny 

Cn  = C nx  + Eny  Dn  = D nx  + Fny 

2.2  Convective  variation 

Let  us  first  recall  that  the  convective  variation  (i.e. 
the  variation  due  to  the  movement  of  the  domain) 


For  a functional  of  the  form  J = / dQ,  the  con- 

vective variation  is: 


where  cj„  is  the  projection  of  w to  the  normal  to 
the  surface  T 

Applying  the  formula  above  to  the  Lagrangian  of 
the  system  (6)  leads  to  the  following  convective 
variation  (where  Gauss’  formula  is  applied) : 

= f V(77  + V’Tfn)-w  + ('H  + V’Tfn)Kdr/ 

Jr, 

+ J [(f  -fv)T  dxip+  (g  -g„)T  dyV1]  undT 

(20) 

If  the  adjoint  equations  hold,  the  variation  of  the 
Lagrangian  is  formed  by  the  convective  variation 
(20)  and  the  local  variation  due  to  the  design  equa- 
tions (8).  Considering  that  Co  = £<5u,  the  optimality 
condition  is: 

sc  = gsu  (21) 

where: 


+ V($  + t£Tfn) 

+ [(f  - f„)T  dx$  + (g  - gt,)T  dyi>]  £n  <ff/] 
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This  equation  is  the  design  equation.  The  quan- 
tity Q can  be  used  as  a gradient  in  an  optimization 
procedure. 

In  summary,  the  optimality  condition  has  the  form 
of  three  sets  of  equations: 

• the  state  equations  (2)  and  boundary  condi- 
tions (3); 

• the  adjoint  equations  (16)  and  boundary  con- 
ditions (17,18,19); 

• the  design  equations  (22). 

3 Application  to  two  dimensio- 
nal Navier-Stokes  equations 


• the  no-slip  condition:  u = v = 0 

• the  isothermal  condition  T — To  or  adiabatic 
condition  VT  • n = 0 


3.2  Optimization  problem 

The  functional  to  minimize  has  the  form: 

J = f $(p,  tw,  rn)dTf  (24) 

Jrf 

with  p the  pressure  and  rw  the  wall  shear  stress. 
The  wall  normal  viscous  stress  rn,  which  is  zero, 
has  to  be  introduced  for  compatibility  reasons  in 
the  adjoint  boundary  conditions  [10] [2]. 


3.3  Adjoint  equations 


The  adjoint  equations  are  given  by  (16).  Expanding 
this  equation  leads  to: 

ATdxij>  + B Tdyi>  - Mt  (dxfv  + dygy^J  = (25) 

with: 


3.1  The  Navier-Stokes  equations 

The  2D  Navier-Stokes  Equations  have  the  form: 

3r(f  - fv)  + 9y(g  - gv)  = o (23) 

where  f,g  are  the  inviscid  fluxes  , and  the  viscous 
fluxes  are  given  by: 

* fv  — [0,  TXxj  Tx y,  /UTXX  T VTXy  hx] 

• [0,  TXy}  Tyy,  UTXy  “f*  VTyy  hy] 

The  components  of  the  stress  tensor  r and  the  ther- 
mal flux  q can  be  expressed  in  term  of  the  first 
derivative  of  the  primitive  variables  q = [p,  u,  v , T\: 

♦ Txx  — %pdxu  — | pdyV 

♦ Tyy  = | fldyV  — | fidXU 

♦ Txy  = ft  ( dyU  + dxv) 

• h = -kVT 


• iv  = [o,  txx , r^y,  kdxfaf 

• gf  — [0,  Tjpy,  Tyy,  kOyl/j^] 

• K - + BjdyVh  or  in  extenso: 

hvi  — 0 

hV2  = TXXdX1p/H  -f"  TXydytj) 4 
hv3  — TXydXll^  4 T Tyy  Oy  1p4 

hV4  = &T  4 • VT^ 

4-  — (Txxdxip 2 + Txydxipz  H-  Txydy4>2  + Tyydyipz) 


where  pt  , are  the  derivatives  of  the  viscosity  and 
conductivity  with  respect  to  the  temperature.  The 
components  of  the  adjoint  tensor  T are: 


with  p(T)  the  dynamic  viscosity,  and  k(T)  the  ther- 
mal conductivity. 

The  far-field  boundary  conditions  are  the  same  as 
for  the  Euler  equations.  The  conditions  for  a wall 


• = ±p  (dx^2  + udxi> 4)  - | {dyil>z  + vdyipi) 

• Tyy  ~ 1 1*  (dy^P  3 + vdyi/>  4)  - § (d*^2  + ufl*^4) 

• T^y  = p (dxi> 3 + dyXp 2 + udy1p4  + VS*^) 


are: 
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3.4  Adjoint  boundary  conditions  on 
a wall 

The  adjoint  boundary  conditions  for  a non-slip  wall 
ar  computed  by  the  equations  (17,18,19).  A lenghty 
calculation  gives  the  following  result: 

• for  ij>: 

~ i’s  = 

- ij) 4 = 0 (isothermal  BC) , or 

kn  ♦ V^4  = $twtw  ^1  — (adiabatic 
BC). 

• for  p: 

- t*n  = + rn) 

_ ,,  — ~ Bx.  _ r 

P$  — Tw  p ->-5 

- fx3  = $twTw  (l  — (isother- 

mal BC) , or 

fi3  = k'ip4  (adiabatic  BC). 

where  ('ipmrfs)  and  (/in,p5)  are  the  normal  and 
t agential  to  boundary  components  of  the  vector 
(*02 j ^3)  and  (jjli,  ^2);  rn  and  Ts  are  the  normal  and 
tangential  components  of  the  adjoint  tensor  T. 

The  adjoint  boundary  conditions  also  give  a com- 
patibility condition  on  the  objective  function 

= -®p  (26) 

3.5  Design  equation 

The  design  equation  is  calculated  using  (22).  Sev- 
eral terms  in  V($  + n)  cancel  each  other  due 
to  adjoint  boundary  conditions,  and  the  gradient  is 
written: 

G = J ^($+p^n  -Twi)s)  dTf 

+ / ($-hp^n  -rw^s)KdTf 

JTf 

+ / pV^n  • 4 - Vi/>s  • £dTf 
Jrf 

+ f pfnVun  * £+ plsVus  • £dV j 
dr  f 

+ Jr  (27) 


with: 

• ^ = i>4TW  (i-f)-  r3 

• /4  = pH'Pi  ~ rn 

4 Discretization  and  numerical 
solution 

The  optimality  condition  is  composed  of  three  equa- 
tions (state, adjoint  and  design).  The  most  com- 
mon method,  which  is  the  one  used  in  this  work, 
is  to  stand  on  the  intersection  of  the  state  and  ad- 
joint surfaces,  i.e.  to  have  state  and  adjoint  vari- 
ables that  satisfy  state  and  adjoint  equations,  and 
to  march  to  the  design  surface  . This  method  is 
the  well-known  descent  method,  commonly  used  in 
mechanical  optimization. 


4.1  Navier-Stokes  equations 

The  Navier-Stokes  equations  (23)  are  discretized 
by  finite- volume  method  and  solved  by  an  implicit 
time  marching  method [9]  . 

Integrating  the  unsteady  equations  on  a finite  vol- 
ume Q i and  using  the  Gauss  theorem  yields: 

dt  f sdQi+  <f  (F(s,n)  — Ft/(s,n))  c?r  (28) 

Jr 

where  F = fnx-\-gny  and  = f vnx  -f  g vny  are  re- 
spectively the  advective  and  viscous  fluxes  through 
the  boundary  T of  the  cell. 

Considering  polygonal  cells,  and  integrating  the  flux 
using  an  n^-points  Gauss  formula,  gives: 

ng 

dt  / S dQi  + 0 E w»  (*i  - **)  = 0 (29) 

Jcii  j£Ni 

where  lj  is  the  length  of  the  edge  jy  wg  the  weight  of 
the  Gauss  point, and  N{  is  the  set  of  edges  defining 
the  control  volume  . The  numerical  inviscid  flux 
are  computed  using  the  Roe  flux-difference  split- 
ting: 

F = I [F(sL,n)  + F(sfi,n)  - \Q(sL,sR,n)\  (s*  -s1-)] 

(30) 
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where  Q is  the  Jacobian  matrix  evaluated  using 
Roe's  average  of  s^s^,  the  latter  being  the  value 
at  the  right  or  the  left  of  the  edges  extrapolated 
from  the  center  of  the  adjacent  cells.  This  recon- 
struction can  be  of  order  0 (constant),  1 (linear) 
or  2 (quadratic),  and  may  imply  complex  detector- 
limiter  to  deal  with  flow  discontinuities  (shocks  or 
slip  lines). 

A centered  discretization  is  used  to  discretize  the 
viscous  terms  of  the  Navier-Stokes  equations.  The 
viscous  flux  is  evaluated  at  one  point  located  at  the 
middle  edge. 

Once  the  spatial  part  is  discretized,  the  equation 
is  integrated  in  time  with  a fully  implicit  Newton- 
like method.  The  following  non-linear  equation  is 
solved  at  each  time  step: 

^(sn+1,sn)  = sn+^~s"  +re(sn+1)  = 0 (31) 

If  this  equation  is  linearized,  the  following  linear 
system  is  obtained: 

(T  + A sn+1  = -^(sVn)  (32) 

with  J = 4^,  the  jacobian  matrix.  This  system 
is  solved  by  a 1-step  Newton  method  using  a ILUO 
preconditioned  GMRES  Algorithm[14]. 

4.2  Adjoint  equations 

The  continuous  adjoint  equations  (16)  are  discretized 
with  the  same  finite-volume  method  as  the  Euler 
equations:  integrating  (16)  over  a cell  i gives  (us- 
ing Gauss  theorem): 

dt  f il>dQi+  f ATdxf  + B Tdyi>dtti 
Jnt  Jcti 

+ f dxf+dygdT=  f h dtli  (33) 
Jfti  Jcii 


where  the  numerical  adjoint  advective  flux  G,j  is 
computed  with  a Roe-like  flux  splitting  formula: 

G.'J  = \ [An(s.)  (V^  + lpR)-l'  (i>R  - ipL)]  (35) 

with  An  — Anr-hBny  and  v an  artificial  dissipation 
factor.  This  dissipation  factor  is  mandatory  for  the 
convergence  of  adjoint  equations  used  in  transsonic 
flows  problems.  It  will  be  shown  that  in  such  prob- 
lems, the  discontinuous  state  variables  imply  Dirac 
peaks  in  the  adjoint  solution.  A lot  of  dissipation 
is  then  needed  to  avoid  divergence  problems. 

As  for  the  state  variables,  the  adjoint  variables  on 
the  edges  are  extrapolated  from  the  nodal  values  us- 
ing a constant,  linear,  or  quadratic  reconstruction. 
The  adjoint  viscous  flux  Fv  is  evaluated  in  a simi- 
lar manner  as  the  viscous  flux  of  the  Navier-Stokes 
equations. 

Time  derivatives  are  added  to  the  adjoint  equa- 
tions to  obtain  a time  dependent  system  of  equa- 
tions, which  is  then  solved  iteratively  by  using  the 
same  method  as  for  the  Navier-Stokes  equations, 
i.e.  fully-implicit  Newton-like  integration. 

4.3  Optimization  procedure 

The  optimization  procedure  presented  here  consists 
in  transforming  the  highly  non-linear  shape  opti- 
mization problem: 

min  /(x) 

with  Ci(x)  < Ci  for  t = 1 • * »nc  (36) 

where  the  function  and  constraints  are  not  explicitely 
known  into  a sequence  of  simpler  (explicit)  prob- 
lems. At  each  iteration  ft,  a subproblem  Pk  : 

min  /(fe)(x) 

with  c^(x)  < C{  for  i — 1 • • - nc  (37) 

is  constructed  and  solved  by  the  CON  LIN  proce- 
dure developed  at  the  University  of  Liege  [6]. 


or,  using  cell  averaged  values  and  applying  Gauss 
theorem  on  integrals  containing  dxip}  dyip, 

ng 

dti>i  + ) ] Ij  y ] Wg  ij  -f  Mj  Fvn  (^ij)! 

j€Ni 

tig 

= A™  £ h £ W9'Pij  (34) 

jeNi 


In  this  procedure,  the  key  point  consists  in  find- 
ing a good  modelling  of  the  objective  function  and 
constraints.  Two  different  approaches  were  tested. 


The  first  method  is  called  ” Method  of  Moving  Asymp- 
totes” (MMA)[llj.  At  iteration  fc,  a function  / is 
approximated  by: 


/<*>(*)  = £ 


— ' rj(k)  — x. 
=1  ui  xi 


■ + ■ 


0* 


- L) 


(*) 


+ r (38) 
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The  parameters  r,p2-  and  g2  are  adjusted  such  that 


/(xW) 

9f  I 


/(x<*>) 

3/ 

dxi 


?(*) 


In  addition  p2  is  set  to  zero  when  dXif  < 0 at 
and  qi  is  set  to  zero  when  dXif  > 0,  such  that  / is  a 
monotoneous  increasing  (resp.  decreasing)  function 
Of 


(. fc ) 

= max 

o, 

(*) 

' = mm 

°,( 

The  parameters  and  L ^ are  the  asymptotes 
which  are  allowed  to  move  from  one  iteration  to 
another,  in  order  to  narrow  or  broaden  the  feasible 
space 


The  second  method  of  approximation  is  the  ’’Di- 
agonal Quadratic”  method  [4].  At  iteration  k,  a 
function  / is  modelled  by: 


PH*)  - /(x(fc))+^3  &,•(*,— *ifc))+^ai  (a;,-  - *jfc)) 

(39) 


with: 


df_ 

dxi  , 

0V 

dxj 


r(*) 


?(*) 


As  the  second  derivatives  of  the  function  is  usually 
not  available,  the  term  a2  has  to  be  approximated. 
We  have  employed  a rank-1  update  procedure  for 
the  coefficients  a2 : 


a^  = a^+aV-L 

‘ * z 


(40) 


with  v = V/(fc+1)  - z = vT  (x(fc+1>  - x<*>) 

and  a a relaxation  factor: 


a = min 


*IMI  ) 


5 Geometry  considerations 

5.1  Airfoil  parametrization 

The  profile  is  defined  by  the  formula: 

y{x)  = yc(x)±yt(x)  (41) 

where  yt  and  yc  are  the  thickness  and  camber  dis- 
tributions defined  as  follow  (see  figures  1(a)  and 
1(b)): 


• the  thickness  distribution  is  modelled  by  two 
Bezier  curves  (one  curve  of  degree  3 on  the 
front  part  and  one  of  degree  2 on  the  aft  part); 
the  design  variables  are  the  position  and  value 
of  maximum  thickness,  the  leading  edge  ra- 
dius and  the  trailing  edge  wedge  angle. 

• the  camber  distribution  is  modelled  by  two 
Bezier  curves  of  order  2;  the  design  variables 
are  the  position  and  value  of  maximum  cam- 
ber, and  the  first  derivative  at  leading  and 
trailing  edges. 


The  airfoil  section  is  then  defined  by  eight  parame- 
ters. 


5.2  Mesh  movement 

During  the  optimization  cycle,  the  geometry  is  mod- 
ified and  the  computational  mesh  must  be  adapted. 
A straightforward  method  consists  in  generating 
a new  mesh  at  each  iteration,  but  this  requires  a 
fast  and  robust  mesh  generation  technique.  In  the 
method  used  here,  which  is  called  ’’spring  analogy”, 
the  vertices  of  the  mesh  are  moved,  without  chang- 
ing the  connectivity  nor  the  number  of  mesh  points. 
The  mesh  is  considered  as  a spring  network  that  re- 
acts to  a perturbation  of  the  boundary  by  adjusting 
the  vertex  position  in  such  a way  that  the  deforma- 
tion energy  is  minimal.  This  leads  to  a system  of 
non-linear  equations  that  are  solved  by  a Newton 
algorithm. 

The  spring  analogy  method  is  well  suited  for  un- 
structured meshes  for  inviscid  flows  but  needs  some 
modifications  to  be  used  on  hybrid  meshes  for  vis- 
cous flows,  because  of  the  high  aspect  ratio  of  the 
cells  in  the  boundary  layer  that  can  produce  over- 
lapping or  highly  distorted  meshes.  A vertex  i near 
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Figure  1:  Airfoil  modelisation 


the  moving  boundary  (in  the  boundary  layer)  is 
moved  according  to  the  formula: 

xr  - x-'d  = a(d)  (x—  - xf  ) (42) 

where  j is  the  nearest  point  on  the  boundary  with 
respect  to  i,  and  a(d)  G [0  1]  is  a decreasing  function 
of  the  distance  to  the  boundary  d,  and  vanishes  for 
d > dcut.  This  procedure  is  applied  to  all  the  ver- 
tices for  which  d < dcut,  and  the  remaining  vertices 
are  updated  with  the  spring  analogy. 


adjoint  solution  The  first  adjoint  variable  obtained 
at  the  first  iteration  step  is  shown  on  figure  2. 


First  adjoint  variable 


(b)  ipi  on  airfoil 


Figure  2:  First  adjoint  variable 

This  peak  is  smoothed  by  the  numerical  computa- 
tion which  involves  a numerical  dissipation  scheme. 
Attempts  to  use  a less  diffusive  flux  (i.e.  a Roe- 
type  flux  difference  splitting)  may  lead  to  a failure 
of  convergence. 


6 Results 

6.1  Inviscid  transsonic  Lift/Drag  Ra- 
tio optimization 

The  test  case  consists  in  maximizing  the  lift  to  drag 
ratio  at  M=0.8  and  an  incidence  of  2 degrees,  with  a 
lower  bound  to  the  maximumthickness  of  the  airfoil 
(12%  of  chord).  In  such  a configuration,  there  is  a 
shock  on  the  lee-side  of  the  airfoil.  This  discontinu- 
ity in  the  state  variable  induces  Dirac  peaks  in  the 


This  type  of  discontinuity  is  inherent  to  the  prob- 
lem, because  the  shock  in  the  flow  solution  is  not 
differentiable.  Such  Dirac  discontinuities  can  be 
found  also  when  using  the  flow  sensitivity  deriva- 
tives to  compute  the  gradient  [8]. 

The  problem  of  minimizing  the  ratio  is  not  obvi- 
ous because  the  cost  function  cannot  be  written  in 
the  form  (24),  although  Ci  and  Cd  have  this  form. 

The  gradient  of  / = ^ is 

V/  = T (vcd  - fvQ) 


(43) 
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One  straightforward  approach  is  to  compute  two 
adjoint  problems,  with  functions  Ci  and  Cd  and 
to  put  their  gradients  in  the  formula  above.  The 
method  used  here  consists  in  using  the  functional 
f ~ ( Cd  ~ foCi),  defined  at  the  current  design 

point  ao,  for  the  computation  of  the  gradient.  This 
function  has  the  same  gradient  as  the  initial  cost 
function  at  the  current  design  point,  and  has  the 
form  (24),  because  it  is  linear  in  Ci  and  Cd-  Its 
gradient  only  requires  one  adjoint  computation. 

After  twenty  iterations,  the  lift  to  drag  ratio  has  in- 
creased from  2.03  * 10-2  to  10.65.  The  final  airfoil  is 
slightly  thinner  than  initial  airfoil,  with  camber  and 
a sharper  trailing  edge,  see  figure  3.  The  lee-side 
shock  has  slightly  increased  and  moved  backward 
(fig.  3). 


(a)  Shape  modification 


(b)  Cp  comparison 

Figure  3:  Shape  and  Cp  comparison 


6.2  Lift  maximization 

The  first  test  case  is  a lift  coefficient  maximization 
of  an  airfoil  in  subsonic  viscous  flow  (Mach =0.5, 
Re=100.000).  The  initial  airfoil  is  a slightly  cam- 
bered NACA0012,  and  an  upper  bound  on  the  drag 
coefficient  is  imposed.  This  problem  was  solved 
with  different  approximations  for  the  objective  and 
the  constraint  ( in  the  text  below,  MMA/QUA  stands 
for  55  Method  of  Moving  asymptots  for  the  objective 
and  diagonal  quadratic  method  for  the  constraint” , 
for  example).  Results  are  presented  in  table  1 and 
2. 


Method 

CL 

CD 

Hsu 
C n 

Initial 

1.778 -lO"1 

3.024  • 10~2 

5.88 

MMA/MMA 

2.399- 10-1 

2.771  • 10-2 

8.66 

MMA/QUA 

2.449  - 10- 1 

2.774  • 10-2 

8.83 

QUA/QUA 

2.273-  IQ-1 

2.945  • 10~2 

7.72 

Table  1:  Aerodynamic  coefficients 


Method 

A CL 

ACd 

MMA/MMA 

+34.9% 

-8.30% 

+47.2% 

MMA/QUA 

+37.7% 

-8.30% 

+50.1% 

QUA/QUA 

+27.8% 

-2.60% 

+31.2% 

Table  2:  Performance  of  the  different 
approximation  techniques 


It  can  be  noticed  that  the  all-MMA  or  MMA/QUA 
strategies  work  better  than  the  all-QUA  strategy. 
That  can  be  due  to  two  facts: 


• the  MMA  approximation  is  well  suited  for 
monotonically  decreasing  or  increasing  func- 
tions (and  that  can  be  the  case  for  Cl); 

• the  QUA  method  relies  on  the  quasi-Newton 
approximation  of  the  Hessian  matrix,  which 
can  be  quite  inaccurate,  especially  during  the 
first  iterations. 


Figure  4 shows  the  initial  and  optimized  shapes  for 
MMA/QUA  and  QUA/QUA  computations  (MMA/MMA 
gives  almost  the  same  airfoil  as  MMA/QUA). 

Figure  5 shows  the  Mach  field  around  the  optimized 
airfoil  (MMA/QUA). 

The  optimized  airfoil  is  cambered  (1.1%)  to  increase 
lift,  with  the  position  of  maximum  camber  displaced 
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Figure  4:  Initial  and  optimized  shapes 


Figure  5:  Mach  field  around  optimized  airfoil  (Mach 
max. =0.89253) 

toward  aft  (at  75%  of  chord),  in  such  a way  that  the 
recirculation  bubble  is  smaller  (see  figure  6),  thus 
leading  to  a decrease  in  the  viscous  drag  coefficient. 


6.3  Inverse  design 

Inverse  design  problems  can  be  treated  as  optimiza- 
tion problems,  the  function  to  optimize  being  the 
least  square  distance  between  a target  pressure  field 
and  the  actual  pressure.  Although  optimization 
techniques  are  much  more  expensive  in  term  of  cal- 
culation time,  they  always  succeed  in  finding  an 
optimum  value  where  inverse  design  techniques  can 
fail.  Inverse  design  problems  are  also  useful  when 
testing  algorithms,  because  the  optimal  solution  is 
known. 


Figure  6:  Skin  friction  coefficient  on  initial  and  op- 
timized airfoils 

In  the  present  test  case,  the  pressure  of  an  air- 
foil (Mach=0.5,  Re=100.000)  has  to  be  matched 
to  the  pressure  field  of  an  RAE2822  airfoil  in  the 
same  flow  conditions.  The  initial  airfoil  is  the  same 
NACA0012  as  in  the  previous  test  case.  The  opti- 
mization technique  used  is  MM  A. 

As  can  be  seen  on  figure  7,  the  final  solution  is  very 
close  to  the  target,  in  term  of  shape  and  pressure 
distribution. 

The  residual  history  gives  valuables  observations  of 
the  behaviour  of  the  MMA  algorithm:  the  initial 
oscillations  of  the  objective  functions,  due  to  oscil- 
lations of  the  variable  describing  camber,  are  pro- 
gressively damped  by  bringing  the  moving  asymp- 
tots  closer  to  each  other.  However  this  leads  to  a 
non-optimal  solution;  this  difficulty  was  overcome 
by  restarting  the  algorithm  after  iteration  13  (see 
figure  8). 

In  the  previous  calculation,  the  upper  and  lower 
bound  on  the  design  variables  are  chosen  such  that 
the  optimal  solution  lies  on  the  boundary  of  the  de- 
sign space.  In  such  a case,  the  MMA  approxiamtion 
is  very  efficient,  because  the  objective  function  is 
generally  a monotonous  function.  When  the  solu- 
tion is  inside  the  design  space,  MMA  fails  to  find 
the  real  optimum  because  it  cannot  model  functions 
with  one  (or  more  ) extremum.  Figure  9 shows 
the  result  of  the  same  test  case,  with  a larger  de- 
sign space:  the  computed  solution  lies  on  the  design 
space  boundary  and  is  not  the  real  optimum. 
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(a)  Shape  comparison 


(b)  Pressure  coefficient  comparison 


Figure  7:  Initial,  target  and  optimized  shape  and 
Cp  comparison 

7 Conclusions  and  prospects 

In  this  paper,  we  have  showed  that  the  adjoint 
method  can  give  accurate  gradient  to  aerodynamic 
coefficients.  This  method,  together  with  fast  flow 
analysis  algorithms  and  a robust  constrained  op- 
timization procedure,  is  mandatory  for  the  resolu- 
tion of  aerodynamic  shape  optimization  problems 
in  2D  on  small  workstations  (all  the  calculations 
were  done  on  a single  CPU  HP-C160  workstation). 

Further  work  is  needed  on  the  resolution  of  the 
adjoint  equations,  especially  to  derive  the  adjoint 
equations  for  turbulent  flows.  Another  topic  of  in- 
terest is  the  study  of  approximations  of  the  func- 
tions; a hybrid  method  involving  both  MMA  and 
QUA  could  be  very  efficient  [12]  if  a good  approx- 


Figure  8:  History  of  objective  function  and  camber 

imation  of  the  Hessian  matrices  could  be  found. 
Further  studies  on  Hessian  matrices[3]  could  thus 
be  valuable. 
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